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Abstract 

Glacial-interglacial cycles are large variations in continental ice mass and 
greenhouse gases, which have dominated climate variability over the Qua¬ 
ternary. The dominant periodicity of the cycles is ^40 kyr before the so- 
called middle Pleistocene transition between ~ 1.2 and ^0.7 Myr ago, and 
it is ^100 kyr after the transition. In this paper, the dynamics of glacial- 
interglacial cycles are investigated using a phase oscillator model forced by 
the time-varying incoming solar radiation (insolation). We analyze the bi¬ 
furcations of the system and show that strange nonchaotic attractors appear 
through nonsmooth saddle-node bifurcations of tori. The bifurcation analy¬ 
sis indicates that mode-locking is likely to occur for the 41 kyr glacial cycles 
but not likely for the 100 kyr glacial cycles. The sequence of mode-locked 
41 kyr cycles is robust to small parameter changes. However, the sequence 
of 100 kyr glacial cycles can be sensitive to parameter changes when the 
system has a strange nonchaotic attractor. 

Keywords: Glacial-interglacial cycles, ice age, quasiperiodically forced 
dynamical systems, strange nonchaotic attractor, SNA, nonsmooth 
saddle-node bifurcations, middle-Pleistocene transition 


1. Introduction 

The glacial-interglacial cycles or simply glacial cycles are alternations 
between cold (glacial) and warm (interglacial) periods, which occurred over 
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the last 3 million years (Myr). They are characterized by large fluctuations 
in the continental ice volume, which can be estimated, among others, from 
the oxygen isotope ratio 5 18 0 of benthic foraminifera sampled in deep-sea 
cores (see Fig. g hhs]. a large value of 5 18 0 indicates a large volume of ice 
sheets. The dominant period of the glacial cycles is ^40 kyr before ^1.2 Myr 
ago and ^100 kyr after ^0.7 Myr ago. This frequency change accompanying 
an increase amplitude is called the middle-Pleistocene transition (MPT) (see 
[6] and references therein). 

There is abundant empirical evidence dHE] that glacial cycles are con¬ 
trolled by the changes in the parameters of Earth’s orbit (eccentricity e and 
longitude of perihelion w) and its obliquity e mi- Milankovitch’s theory pro¬ 
vides a physical mechanism for this control |12| : the orbital and obliquity 
changes determine the seasonal and spatial distributions of insolation (this 
is called the astronomical forcing ), and, specifically, northern hemisphere 
summer insolation controls the interannual accumulation of snow and the 
subsequent growth of ice sheets. Milankovitch’s mechanism is nowadays 
viewed as realistic, but it is considered to be only one element of a com¬ 
plex nonlinear dynamical process. For example, the glacial cycles of the last 
0.7 Myr have a dominant spectral power near 100 kyr, while the astronom¬ 
ical forcing has negligible power near 100 kyr (so called the 100-kyr cycle 
problem 13]). Thus, some nonlinear transformation mechanism must exist 
from astronomical forcing to glacial cycles [7j. 

The glacial cycles have been investigated using various models rang¬ 
ing from simple conceptual models to complex coupled general circulation 
models (for example, |13|). between which there are Earth Systems Models 
of Intermediate Complexity (EMICs) [HHT6]. Among others, conceptual 
models provide theoretical frameworks to identify the dynamics of glacial 
cycles (for example, see [J, T7P29] and references therein). They are pre¬ 
sented under the form of low-dimensional dynamical systems forced by a 
measure of northern hemisphere insolation or, more generally, by a linear 
combination of esinzu (the climatic precession ) and obliquity e. Since these 
quantities are well-approximated as quasiperiodic functions of time over the 
past several million years mi, the models of glacial cycles may be viewed 
as quasiperiodically forced dynamical systems. 

Let T n = R n /(2ttZ) n be an TV-dimensional torus. Quasiperiodically 
forced dynamical systems can be represented in a skew-product form: 

(e = u, 6 e t n , 

1 x = f(x, 9), xEl M , 
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where 6 — (0i, 62 ,0w) T is the phase of the drive system, x = (aq, aq,^ m) T 
is the state of the response system, f (x, 6) is a periodic function in each phase 
Oi , and uj — (uji,uj 2 , ...,c^w) T is a vector of incommensurate frequencies such 
that feiCJi + ^ 2^2 + • • • + = 0 does not hold for any set of integers, £q, 

& 2 , fc/v, except for the trivial solution £q = /c 2 = *** = &/v = 0. In the 
models of glacial cycles, 0(£) corresponds to the phase of the astronomical 
forcing, and x(£) corresponds to the climate state. 

Quasiperiodically forced systems can exhibit intermediate dynamics be¬ 
tween quasiperiodicity and chaos, so-called strange nonchaotic attractors 
(SNAs) 32. 33] . An SNA is a geometrically strange attractor for which typ¬ 
ical Lyapunov exponents are nonpositive [32j (see pa Ea for comprehensive 
reviews). The system 0 is characterized by (TV + M) Lyapunov exponents. 
Among them, N Lyapunov exponents are trivially zero. They correspond to 
the phase equations of the drive system. When the system has an SNA, the 
nontrivial largest Lyapunov exponent is negative in general. Thus, under a 
common quasiperiodic forcing, trajectories x(£), which start from different 
initial conditions x(0), approach (or synchronize to) a unique, or one of a 
finite number of possible trajectories as time elapses (this is a synchronizing 
property) [3lJ[35]. However, related to the strange geometry of SNAs, tra¬ 
jectories x(£) have nonexponential sensitivity on initial phases 0(0) [361 139] 
or on parameter values of the systems m- SNAs have been observed in 
many laboratory experiments [ 35 , EHH)] but very rarely in nature so far. 
Linder et al. show that the brightness changes of some RRc Lyrae stars 
have nonchaotic and fractal properties of SNAs m ■ 

Recently, the authors of this paper showed that several models of glacial 
cycles exhibit SNAs [23 EH]. This means that the relationship from the 
phase of astronomical forcing 6 to the climate state x is represented by a ge¬ 
ometrically strange set. When a model of glacial cycles exhibits an SNA, the 
sequences of glacial cycles synchronize under the same astronomical forcing, 
but simultaneously they can be sensitive to parameter changes. Ivashchenko 
et al. showed that the autocorrelation function of the benthic isotopic data 
has self-simility characteristic of SNAs [48] . 

However, so far, the birth mechanism of SNAs has not been elucidated 
in the models of glacial cycles. In this paper, we introduce a phase oscillator 
model of glacial cycles, whose bifurcation analysis is easier than higher¬ 
dimensional models, and show the birth mechanism of SNAs in this model. 
Based on the bifurcation analysis of the phase model, we suggest that the 
41-kyr cycles, typical before the MPT, are little sensitive to the parameters. 

By contrast the sequence of 100-kyr cycles, which characterises climate after 
the MPT, can be strongly sensitive to small parameter changes. 
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The remainder of this article is organized as follows. In Section 2, a phase 
oscillator model for glacial cycles is introduced. In Section 3, the bifurcations 
exhibited by the model are analyzed for two types of external forcing: an 
ideal two-frequency quasiperiodic forcing and the astronomical forcing. In 
Section 4, we discuss parameter sensitivity of the phase model, and the 
MPT. Section 5 summarizes this article. In Appendix, we characterize the 
strangeness of attractors using the phase sensitivity exponent [36] . 

2. Model 

Consistently with the above discussion, the astronomical forcing F^it) 
is calculated by the following formula [26] : 



where we set a = 23.58 W/m 2 so that Ta(£) has unit variance. The fre¬ 
quencies uoi and coefficients, and q, are based on [26], and they are listed 
in m in descending order of power sj + c 2 . The first three periods corre¬ 
spond to the climatic precession: 2tt/oji ~ 23.7 kyr, 2n/uj 2 ~ 22.4 kyr, and 
27r/cc;3 « 19.0 kyr. The fourth period 2n « 41.0 kyr corresponds to the 
obliquity change m- Unless mentioned otherwise, we use the time unit of 
10 kyr and the angular unit of radian for the variables and parameters of 
the model. 

Glacial cycles are described here by a forced oscillator, as suggested by 
some recent analyses mm- Further, the global ice volume V (0) is modeled 
as a 27r-periodic function of the phase of the oscillator (j){t) E T 1 = R 1 / (27tZ), 
i.e., V((f>) — V((f> + 2tt). We express the ice volume U(</>) as a Fourier sine 
series up to the second harmonic: 

V ((f>) =—sin (f) —-sin 2(f), (3) 

where S is the modification parameter. For |<5| < 1, the ice volume U(</>) has 
a unique minimum at </> m j n = cos _1 [(v / l + 8 S 2 — 1)/(4<5)] and a maximum at 
(f> max = 2vr - (f) min , as shown in Fig. (2^a). 

In general, the phase motion of a limit-cycle oscillator under a weak 
forcing p(t) is described by </> = /? + Z(0) • p (£), where f3 is the natural 
frequency of the oscillator, and Z(</>) is the sensitivity function with Z(</>) = 
Z(0 + 2n) [49] . Based on this general phase equation, we suppose the phase 
motion of glacial cycles as follows: 

j> = p-aV'(<t>)[l + 'yF(t)], (4) 
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where U'(0) — dV }^ is the gradient of the ice volume, F(t) is a forcing 
function such as i%(£), and the parameters cq 0, and 7 are assumed to be 
positive. The specific term, —aV'(</>)[ 1 + 7 F(t)], provides a forcing to move 
toward an ice volume minimum for a large astronomical forcing with F(t) > 
— 1/7 and to move toward an ice volume maximum for a small astronomical 
forcing with F(t) < — 1 / 7 . We can also write Eq. (4) as 0 = — dU Q^ using 
a time-dependent potential 1 /( 0 %) = — 00 + aV(q>) 1 + 7 E(£)] (Fig. [ 2 ^a)). 
In this paper, we restrict ourselves to the simple case \S\ < 1/4, where the 
number of local minima of 17(0, t) as a function of 0 is at most one. 

Let us first consider the case of constant forcing F(t) = F c is considered. 
The bifurcation diagram of Eq. Q for F(t) = F c is shown in Fig^b). For 

^, a stable equilibrium point, or a stable 
, and an unstable equilibrium 
Thus, the system approaches 


low insolation F r < — 


if- 

7 V c 


+ 1 


^a(l+<5) 

node , exists near the ice volume maximum 
point, or a saddle , near the minimum 


there is 


a stable glacial state. For -I ( a J j +S) + l) < F c < I f a( /7) “ 1 )> 
no equilibrium point, and glacial-interglacial oscillations emerge. For high 


insolation F c 


>i 


0(1-5) t)’ 


a(l-S) 

and a saddle near the maximum 


a stable node exists near the ice volume 
ax. Thus, the system 


minimum <p n „„, 
approaches a stable interglacial state. 

Next, the case of astronomical forcing F(t) — J%(t) is studied. Unless 
otherwise noted, differential equations are integrated by the Runge-Kutta 
4th order method with a step size of 100 yr (namely h — 0.01). Figure [5] 
shows a simulated ice volume trajectory V ( 0 (t)) for a = 7 = 1 . 0 , 0 = 1.0006, 
and 5 = 0.24 with initial condition 0(£o) = 0 at % = —20 Myr. Here, 
assuming a = 7 = 1.0 and 6 = 0.24, the value of 0 was tuned to maximize 
the Pearson’s correlation coefficient r over the past 700 kyr between the 
ice volume reconstruction [ 5 ] and a simulated ice volume trajectory U( 0 (£)) 
with 0(£o) = 0 at to = —20 Myr. This yields r — 0.76 for 0 = 1.0006. 

Generally, trajectories 0(t) can depend on initial conditions 0(to)* How¬ 
ever, due to the astronomical forcing i%(t), trajectories starting from differ¬ 
ent initial phases 0 (to) can converge to a single or some pieces of trajectories 
after a certain transient time, as shown in Figs.^a) andQb). The transient 
time can last for several million years for the case of 100 -kyr cycles, but it 
is only several cycles for the case of 41-kyr cycles. 

It should be noted that the asymmetric sawtooth pattern of the ice 
volume trajectory U( 0 (£)) is not due to the subtle asymmetry of function 
U(0) caused by the second harmonic —(S/2) sin20 but due to the specific 
form of Eq. ©; the phase 4>{t) increases slowly in glaciation phases and 
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quickly in deglaciation phases, on average. The qualitative features of the 
following bifurcation analysis do not change even if the second harmonic is 
omitted or changed in the range \S\ < 1/4, but we keep the term because it 
increases the agreement with data, for example, by about 10% from r = 0.64 
for S — 0 to r — 0.76 for 5 — 0.24. 

3. Bifurcation analysis 

We study bifurcations of the model for the astronomical forcing Ta(£) 
and an ideal two-frequency quasiperiodic forcing described below. We set 
7 = 1.0 and S = 0.24, and regard a and /3 as control parameters. 

3.1. The case of the ideal two-frequency quasiperiodic forcing 
We start from the two-frequency quasiperiodic forcing 

Fc(t) = sincjit + sine(5) 

where uo\ is the most prominent frequency of the climatic precession ( 2 tt/uji ~ 
23.7 kyr), and ujg = x (y/E — l)/2 (2 tt/ujg ~ 38.4 kyr) is a frequency near 
the dominant frequency of obliquity change 2 ^k joj 4 ~ 41.0 kyr nu, which 
was chosen to be most incommensurate with uj\. Defining phase variables, 
61 = uj\t (mod 2 i r) and 6 q = co G t (mod 2 i r), the model can be expressed in 
the skew-product form 

fj) = ^ + (a(cos0 + 5cos20)[l +jF g ( 0 i, 9 g )\, 

( 6 ) 

= CUI, 0g = UG, 

where F G {t) is redefined as F G (0 1 , 6 G ) = sin^i + sin 9g, and the phase 
space is T 3 . Similar pendulum-type systems but with additive quasiperiodic 
forcing have been studied in j50H53] . The attractor of the system is obtained 
by solving Eq. ^ forward for a long time. In the particular case of this 
model, the unstable invariant set of the system is found by solving Eq. © 
backward for a long time (this is called the repeller) [52]. We plot the 
attractor and the repeller at a Poincare surface of section (PS) at a certain 
constant value of #i(~ 2.168184) Q 

The system ^ has three Lyapunov exponents: one is the Lyapunov 
exponent in the direction of 0, and the other two are trivially zero, 


Tn numerical calculations for Fig. [H] a step size h = 2n/uj\ x 0.002 was used so that 
6 1 takes a constant phase in every 500 steps. 
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corresponding to 6 \ — uj\ and Qq — uq- The dynamics of Eq. ^ are 
classified by the Lyapunov exponent A^ and the winding number W. The 
winding number W is defined as 

W = lim 

t—yoo t 

for the same Eq. ©> but with <p on the line R. The dynamics are referred to 
as mode-locked if the winding number W is constant with respect to slight 
changes of /3 [54]. When the dynamics are mode-locked, the model has a 
two-frequency quasiperiodic attractor and a two-frequency quasiperiodic re- 
peller, which are represented by a stable and an unstable invariant curve 
in PS, respectively (Fig. [5^a) for example), and the winding number W is 
rationally related to the forcing frequencies, i.e., W — (k/m)uji + ( l/m)uG 
(k,l G Z, m G N) [54]. The denominator m corresponds to the multiplicity of 
each invariant curve in the direction of <j>m- The two-frequency quasiperi¬ 
odic attractors are characterized by the negative Lyapunov exponent < 0. 
When the dynamics are non-mode-locked, we find three-frequency quasiperi¬ 
odic motion when A^ = 0 and SNAs when A^ < 0, as is the case of the general 
quasiperiodically forced circle map m (see Figs. m andj^f)). 

The regime diagram for Fo(t ) is shown in Fig. [6^a), which is obtained 
by taking a grid size of A a = 0.01 and A/3 = 2.5 x 10 -3 . The dynamics are 
numerically regarded as mode-locked if the change of W is less than 10 -4 
when (3 is changed either by ±A/3 (white regions). The non-mode-locked 
regions are divided into regions with |A^| < 5.0 x 10 -4 (green) and regions 
with A^ < —5.0 x 10 -4 (magenta). Under these numerical criteria, we find 
typically the two-frequency quasiperiodic attractors in mode-locked regions 
(white), three-frequency quasiperiodic motion in non-mode-locked regions 
with |A^| < 5.0 x 10 -4 (green), and SNAs in non-mode-locked regions with 
< —5.0 x 10 -4 (magenta). The boundaries of each regime depend on the 
grid size and thresholds slightly^] but the overall structure is fairly robust. 
The strangeness of attractor in the non-mode-locked regions with neg- 


2 For the (N + l)-frequency quasiperiodic motion with the Lyapunov exponent of 

zero, the order of the computed Lyapunov exponent converges to zero inversely propor¬ 
tionally to the averaging time T of the local Lyapunov exponent. Typically, we have 

|A</,| ~ O(10 -5 ) for T — 10 2 * * 5 . Thus, the threshold at A^ = —5.0 x 10 -4 is effective to 
distinguish the regimes with and without the Lyapunov exponent A^ of zero. The changes 

in the threshold between —10 -5 < A^ < —10 -4 shift the boundaries between (TV + In¬ 
frequency quasiperiodic motion (green) and SNA (magenta) in the direction parallel to 
nearby mode-locked regions. This uncertainty of the boundaries is typically less than 0.1 
with respect to a and j3. The choice of the threshold for the winding number W is also 
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ative Lyapunov exponent (magenta) may also be assessed with the phase 
sensitivity exponent [36] . which characterizes the sensitivity of with re¬ 
spect to the changes in 0\ or 9q. The phase sensitivity exponent is computed 
in the Appendix, and the results show that the attractor is smooth in the 
mode-locked regions, and strange in the non-mode-locked regions with neg¬ 
ative Lyapunov exponent. 

The mode-locked regions form Arnol’d tongue-like structures in the pa¬ 
rameter plane, but the width of each mode-locked region first increases and 
then decreases to values close to zero. Such leaf-like Arnol’d tongues are 
characteristic of mode-locking in quasiperiodically forced systems m , and 
they have been found in several models of glacial cycles mm- 

The system exhibits two types of bifurcations on the boundaries of mode- 
locked regions. The smooth saddle-node bifurcation of tori occurs on the 
boundaries where the Lyapunov exponent A^ changes from negative to zero 
(between white and green regions). The nonsmooth saddle-node bifurcation 
of tori m F55l |56] occurs on the boundaries where A^ remains negative 
(between white and magenta regions). Here, “tori” means the quasiperiodic 
attractor and repeller in the mode-locked regions (In higher-dimensional 
systems, the latter need not be a repeller but an unstable quasiperiodic 
torus of saddle type). In this paper, we focus on the bifurcations for the 
mode-locked region with W = 0, but they occur in the other mode-locked 
regions as well in the same manner. 

In the smooth saddle-node bifurcation, the distance between stable and 
unstable invariant curves decreases to zero for every value of 0 g, as shown in 
Fig- 05b). At the bifurcation point (3 = /3 C , the stable and unstable invariant 
curves merge at every point in Oq yielding three-frequency quasiperiodic 
motion, as shown in Fig. [5](c) . Specifically, define the distance g?(0i,0g) 
between a point on attractor and a point on the repeller at the same ( 0 i, 6 q), 
and consider its maximum d max and its minimum d m i n on (6 Mg) GT 2 [57]. 
As shown in Fig. [7^ a), they obey a scaling law d m - m C\ \/3 c — (3 1 0 ' 5 and 
dmax — C2\/3 c — /3| 0 ' 5 (Cl < C 2 ). This scaling clearly holds with Ci = C 2 in 
the small forcing limit 7 -T 0 , where the smooth saddle-node bifurcation of 
tori degenerates into that of equilibrium points. 

On the other hand, in the non-smooth saddle-node bifurcation, the colli¬ 
sion of the stable invariant curves and unstable invariant curves occurs only 
in a countable dense set of Oq, as demonstrated in Fig. [5](e) , and it creates 


based on the numerical convergence speed. The above choice gives consistent results with 
the analysis of the phase sensitivity exponent shown in Appendix. 



an SNA and a strange repeller. An SNA and a strange repeller for a — 2.2 
and /? = 1.4 are shown in Fig. [5^f). At the bifurcation point /? = /3 C , the 
minimal distance d m i n decreases to zero, but the maximal distance d max is 
strictly positive, as shown in Fig. Qb) m\ ■ In the nonsmooth saddle-node 
bifurcation, the distances do not show the scaling low \/3 c — /3| 0 ' 5 , which 
appears in the smooth saddle-node bifurcation. 

3.2. The case of the astronomical forcing 

We now consider the case of the astronomical forcing FA(t), where Eq. ^ 
is defined using Fa( 0 i, # 35 ) = ^ sin 6 ^ + q cos^), Qi — c instead 

of Fq(0 1, 6q )• The regime diagram for FA(t) is shown in Fig. [6^b), which is 
produced in the same manner as Fig.^a). As a natural extension from two- 
frequency forcing, the system exhibits a 35-frequency quasiperiodic attractor 
in each mode-locked region (white), 36-frequency quasiperiodic motion in 
non-mode-locked regions with = 0 (green), and an SNA in non-mode- 
locked regions with A^ < 0 (magenta). The strangeness of attractor in 
non-mode-locked regions with negative Lyapunov exponent is assessed by 
the phase sensitivity exponent m in Appendix, which numerically shows 
that the attractor is smooth in the mode-locked regions and strange in the 
non-mode-locked regions with negative Lyapunov exponent. 

As is the case of two-frequency forcing, smooth and nonsmooth saddle- 
node bifurcations occur at the boundaries of mode-locked regions. Due 
to high dimensionality, the attractor and repeller are not visible but we 
can calculate the maximal distance d max and the minimal distance d m i n 
between the attractor and the repeller. On the boundaries with zero Lya¬ 
punov exponent, the smooth saddle-node bifurcation occurs creating 36- 
frequency quasiperiodic motion, with the scaling law d m i n ~ Cl\Pc-P I 0 - 5 
and dmax — C2|/3c — /?| 0 ' 5 , as shown in Fig. (7^c). On the boundaries with 
negative Lyapunov exponent, the nonsmooth saddle-node bifurcation occurs 
creating an SNA and a strange repeller. At the bifurcation point (5 — /3 C , 
the minimal distance d m i n decreases to zero, but the maximal distance d max 
is strictly positive, as shown in Fig. [T](d) . In the nonsmooth saddle-node 
bifurcation, the distances do not show the scaling low \/3 c — /3| 0 ' 5 , which 
appears in the smooth saddle-node bifurcation. 

4. Discussions 

4.1. Parameter sensitivity 

A number of previous studies report parameter sensitivity of the se¬ 
quence of glacial cycles generated by simulation models mmmmmr 
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1291 58]. In some models, parameter sensitivity is attributed to thresholds 
with discontinuity p], and in other models, it is attributed to a chaotic prop¬ 
erty, that is, the exponential sensitivity to initial conditions [29] . However, 
parameter sensitivity can appear even when the models are smooth dynam¬ 
ical systems and nonchaotic. An example is Saltzman-Maasch (1990) model 
[59] as pointed out in |28j. Our phase oscillator model also shows such a 
parameter sensitivity. 

To indicate “the region of 100 kyr cycles” in parameter space, we again 
calculate the Pearson’s correlation coefficient r over the past 700 kyr between 
the ice volume reconstruction [5] and a simulated ice volume trajectory with 
</>(£o) = 0 at to = —20 Myr. The parameter points with high correlation 
r > 0.7 are plotted by symbol “+” in Fig. |6^b)j^] They are distributed in 
a narrow region, where the winding number W is constrained between 0.58 
and 0.66. Around the region of 100 kyr cycles, the mode-locked states are 
unlikely to occur because they have little measure. Even if the mode-locked 
states occurred, they would be fragile against small parameter changes. 

If the 100 kyr cycles are non-mode-locked, the dynamics of SNA or the 
dynamics of (N + 1)-frequency quasiperiodic motion with the Lyapunov ex¬ 
ponent of zero are candidates for the 100 kyr cycles. The dynamics of SNAs 
show the synchronizing property under the same astronomical forcing, as 
shown in Fig. Qa), which is often considered as an important feature for 
the models of glacial cycles. However, in the SNA regime, the sequence of 
glacial cycles can be highly sensitive to parameter changes. Two simulated 
ice volume trajectories for two slightly different parameter values /? = 1.0006 
(red) and (3 = 1.0001 (blue) are shown in Fig. [8](a) . The deviation of the 
trajectories with nearby parameters is reminiscent of the well-known but¬ 
terfly effect in chaotic systems but this system is nonchaotic and even has a 
negative Lyapunov exponent for these parameter values of /?. Since SNAs 
have strange geometrical structures as shown in Fig. [5^f), subtle changes of 
parameter values can induce large shifts of the points in the attractors. 

On the other hand, the 41 kyr cycles is caused by the mode-locking to 
the 41 kyr component of astronomical forcing (due to obliquity motion) in 
this model. The mode-locking region of the 41 kyr cycles is relatively wide 
(see the white region labeled by W — 004 in Fig. [6^b)). Thus, this mode¬ 
locking is likely to occur. This is consistent with the previous study m , 


3 High values of correlation coefficient r with r > 0.7 can be obtained in regions with 
36-frequency quasiperiodic motion (green) or in mode-locked regions with m > 2 (white) 
by chance because trajectories depend on initial conditions over the past 700 kyr. 
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which proposes that the phase of 41 kyr cycles is locked by the astronomical 
forcing. 

The sequence of the mode-locked 41 kyr cycles is relatively robust to 
parameter changes. Figure^b) shows ice volume trajectories corresponding 
to different values of f3 in the mode-locked region: (3 — 1.61, (3 = 1.65, and 
/3 = 1.69. The other parameters are a = 0.8, 7=1, and S — 0.24. 

4-2. Middle Pleistocene transition 

The middle Pleistocene transition (MPT) began ^1.2 Myr ago and was 
complete by ^0.7 Myr ago, through which the average period of glacial 
cycles changed from ^40 kyr to ^100 kyr, accompanying an increase of 
amplitude | 6 ]. A full investigation of the causes and dynamics of the MPT 
is beyond the scope of this paper. In particular, a phase model provides 
no information on amplitude. Yet, we note that a frequency change similar 
to the one observed during the MPT is obtained when the parameter f3 is 
changed over the period of MPT, as shown in Fig. [TJ The parameter (3 is 
fixed at f3 — 1.65 in the model-locked region with W — 004 until 1300 kyr 
ago (see Fig. | 6 ^b)). Then, (3 is decreased linearly until it reaches f3 = 0.9 at 
600 kyr ago. After 600 kyr ago, f3 is kept at [3 = 0.9. The other parameters 
are set at a — 0.8, 7 = 1.0, and S — 0.24. The sequence of the glacial cycles 
in Fig. [l] is robust to changes in initial conditions because of the strong 
synchronizing property of the 41-kyr cycles before the MPT (cf. Fig. |4|(b)) . 

5. Summary 

We introduced a phase oscillator model of glacial cycles and analyzed the 
bifurcations of the model for the ideal two-frequency quasiperiodic forcing 
and for the astronomical forcing. It was shown that SNAs appear through 
nonsmooth saddle-node bifurcations of tori in the model. Based on the 
results for the phase oscillator model, we conjecture that the bifurcations 
from quasiperiodic attractors to SNAs found in oscillator models of glacial 
cycles [28] are also nonsmooth saddle-node bifurcations. The regime diagram 
in Fig. [ 6 ](b) indicates that mode-locking is likely to occur for the 41 kyr 
glacial cycles but not likely for the 100 kyr glacial cycles. The sequence of 
mode-locked 41 kyr cycles is robust to small parameter changes. However, 
the sequence of 100 kyr glacial cycles can be sensitive to parameter changes 
when the system has an SNA. 

Long transient dynamics of million-year scale can be observed for the 
100-kyr glacial cycles though it can vanish quickly for the 41-kyr cycles 
(Figs. m and Q b)). Given that the Quaternary ice age is ~3 Myr, and 
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that there can be perturbations to the system, transient trajectories may be 
important to understand the glacial cycles, not only just on the attractor. 
This problem will be explored elsewhere ED]. 
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Appendix. Strangeness of attractors in the non-mode-locked re¬ 
gions with the negative Lyapunov exponent 

We present numerical support for the existence of SNAs in the non¬ 
mode-locked regions with the negative Lyapunov exponent < 0 (magenta 
regions of Figs. [6^ a) and (b)). Denote formally the phase equation of an 
oscillator under an A-frequency forcing by 


0 = cj, 

4 > = /(<M), 


(7) 


where 9 = (# 1 , 62 , 0 jv) T and u = (car, u> 2 , cnjv) T . When the system Q 
has an A-frequency quasiperiodic attractor, it is represented by a certain 
smooth function <fi — H{6). On the other hand, when the system Q has an 
SNA, the relationship between 6 and </> is discontinuous. 

The phase sensitivity function T{t) of variable (j) with respect to variable 
0 i is defined as 


T(t) = min 

MO), 0(0)} 


< max 

[0<T<t 



( 8 ) 


This derivative 90(r)/90j is obtained by integrating 


d_ fdf \ = 9) fdf\ df( 0, 9) 

dt \d9i) dct> \d9i) d9 % 


(9) 


along the solution of Eq. ©• When the system 0 has an A-frequency 
quasiperiodic attractor, the state (#(£), </>(£)) and the quantity d(j){t)/d0i, 
which started from some initial conditions, approach the attractor and its 
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derivative H f (Q(t )), respectively, as time t increases [36J. Thus, T{t) is 
bounded for ^-frequency quasiperiodic attractors. If we write the phase 
sensitivity function formally as T{t) ~ t^ for large £, the exponent /i is zero 
for ^-frequency quasiperiodic attractors. The exponent [i is called the phase 
sensitivity exponent. On the other hand, SNAs do not have finite derivatives 
by definition, but we can calculate the quantity d<j){t)/dQi (and T(t)) along 
orbits approaching an SNA. It is known that the phase sensitivity function 
grows as T(t) ~ with fi > 0 for SNAs (see 3T. 36J for detail). 

For the case of two-frequency forcing Fq{Q\, Qq), we calculate the phase 
sensitivity function T(t) with respect to Oq , changing parameter /?, where 
the other parameters are set to a = 2.2, 7 = 1, and S = 0.24, as in Fig. [T](b) . 
For the case of astronomical forcing Fa(0\, ..., # 35 ), the phase sensitivity 
function T(t) is calculated with respect to # 4 , changing parameter /?, where 
the other parameters are set to a = 1.4, 7=1, and 5 — 0.24, as in Fig. [7](d) . 
The phase sensitivity exponent /a is measured from the growth of the phase 
sensitivity function over the time interval [ 10 5 , 10 7 ]. 

Figures [9^ a) andj^b) show the winding number W, the phase sensitiv¬ 
ity function ju, and the value of phase sensitivity function T at t — 10 7 as 
functions of f3. Mode-locked regions and non mode-locked regions are shown 
by thick points and thin points, respectively, in the graph of the winding 
number W. The phase sensitivity exponent fi takes a positive value enough 
away from zero in the non-mode-locked regions. This means that the at¬ 
tractors in the non-mode-locked regions are typically strange. On the other 
hand, the phase sensitivity exponent [i takes the value of zero in most of the 
mode-locked regions, which indicates Af-frequency quasiperiodic attractors. 
In some regions (1.4075 < /? < 1.41 and 1.4225 < /? < 1.435 in Fig. ^a)), 
the numerical value of /i obtained with a simulation length of 10 7 can be 
positive due to slow convergence. However, the value of the phase sensitivity 
function T( 10 7 ) itself is small when compared to non-mode-locked regions. 
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Figure 1: (Color) Glacial cycles over the Quaternary: The global ice volume reconstructed 
from the oxygen isotope ratio (black solid line) [5] and the simulated ice volume V((j)(t)) 

(red dashed line). The former is scaled as 1.3 x S ls O — 4 for comparison. The parameter 
ft is fixed at ft = 1.65 until —1300 kyr and fixed at (3 — 0.9 after —600 kyr, between which 
P is decreased linearly interpolating the two values. The other parameters are a = 0.8, 

7 = 1.0, and 5 = 0.24. The sequence is robust to changes in initial conditions because 
of the strong synchronizing property of the 41-kyr cycles before the middle-Pleistocene 
transition. 
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Figure 2 : (Color) (a) Ice volume V((j>) (black solid line) and potential U((j),t) as functions 
of 4> for F(t) =3 (red solid line), F(t) =0 (magenta dashed line), and F(t) == —3 (blue 
dotted line), a = 1.0, f3 — 1.0006, 7 = 1.0, and 5 — 0.24. (b) Bifurcation diagram of 
Eq. @ for constant forcing F(t) = F c . The positions of stable nodes and saddles are 
presented by black solid lines and green dotted lines, respectively. The horizontal lines 
are drawn at the ice volume maximum 0 — 0 ma x and minimum 0 — 0 m in. 



Time [kyr] 


Figure 3: (Color) Global ice volume reconstructed from benthic oxygen isotopic ratio 
(dashed line) [5] and the simulated ice volume V (blue line) using Eqs. © and Q. 
a = 1.0, /3 = 1.0006, 7 = 1.0, and 5 = 0.24. Both curves are normalized. 
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Figure 4: (Color) Synchronization of ice volume trajectories V(cf)(t)) started from different 
initial conditions <j>(to) = 2irk/b (k = 0, 1, 2, 3, 4) at time to = —10 Myr under the 
astronomical forcing Fa (t ). (a) The case of 100-kyr cycles for a = 1.0, = 1.0006, 

7 — 1.0, and 5 — 0.24. (b) The case of 41-kyr cycles for a = 0.80, f3 m 1.65, j — 1.0, and 
6 = 0.24. 
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Figure 5: (Color) Poincare surface of section plots at a certain constant value of 0 1 (« 
2.168184) obtained by integrating Eq. ^ forward (red) and backward (blue) for a long 
time, respectively. Each panel corresponds to different values of parameters a and /3: 
(a) in two-frequency quasiperiodic regime, (b) same but just before the smooth saddle- 
node bifurcation occurring at a — 0.4, /3 C ~ 0.303754, (c) in three-frequency quasiperiodic 
regime, (d) in two-frequency quasiperiodic regime, (e) same but just before the nonsmooth 
saddle-node bifurcation occurring at a = 2.2, f3 c « 1.3964061, and (f) in the SNA regime. 
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Figure 6: (Color) (a) Regime diagram for the two-frequency forcing Fc(t) (N = 2) and 
(b) regime diagram for the astronomical forcing Fa($) (N — 35): the mode-locked regions 
corresponding to iV-frequency quasiperiodic attractors (invariant curves in PS) (white), 
the regions with |A</>| < 5.0 X 10 -4 typically corresponding to (N + l)-frequency quasiperi¬ 
odic motion (green), and the regions with < —5.0 X 10 -4 typically corresponding to 
SNAs (magenta). Prominent mode-locked regions are shown with their rotation numbers. 
The parameter points corresponding to correlation coefficient r larger than 0.7 are marked 
by symbol “+” (see text for the details of the calculation of r ). 
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Figure 7: (color) The minimal distance d m i n and the maximal distance d ma x between the 
attractor and the repeller as functions of /3: (a) a smooth saddle-node bifurcation along 
a = 0.4 for the two-frequency forcing Fc(t). Inset is a log-log plot showing the scaling 
oc | p c — /3\ 0 ’ 5 just before the bifurcation occurring at f3 — /3 C . (b) a nonsmooth saddle- 
node bifurcation along a = 2.2 for Fc(t). (c) a smooth saddle-node bifurcation along 
a = 0.4 for the astronomical forcing Fa(£). (d) a nonsmooth saddle-node bifurcation 
along a = 1.4 for Fa(£). The panels (a) and (c) show the results only for the region with 
attractor (3 < /3 C . 



Time [kyr] 

Figure 8: (Color) (a) Simulated ice volume trajectories for two slightly different values 
in SNA regime: f3 — 1.0006 (red solid line) and = 1.0001 (blue dashed line) (a = 1.0, 
7=1, and S = 0.24). (b) Simulated ice volume trajectories for three different values 
in 41-kyr mode-locked region : f3 — 1.61 (red solid line), f3 = 1.65 (green dashed), and 
f3 — 1.69 (blue dotted line) (a = 0.8, 7=1, and S = 0.24). 
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Figure 9: The winding number W, the phase sensitivity exponent /x, and the value of 
the phase sensitivity function T at t = 10 7 as functions of parameter /3, respectively, 
(a) The case for the two-frequency forcing Fc(t) with parameters a = 2.2, 7 = 1, and 
5 = 0.24. The phase sensitivity function is calculated with respect to 0g • (b) The case for 
the astronomical forcing Fa(F) for with parameters a = 1.0, 7 = 1, and S = 0.24. Mode- 
locked regions and non mode-locked regions are shown by thick points and thin points, 
respectively, in the graph of the winding number W. The phase sensitivity function is 
calculated with respect to # 4 . The exponent /x is calculated from the growth of the phase 
sensitivity function over the time interval [10 5 , 10 7 ]. The phase sensitivity function is 
minimized for 10 initial conditions. 
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